;+
; NAME:
;       BGPS_POSTAGE
;
; PURPOSE:
;       Creates postage-stamp images of each BGPS source, centered on
;       the pixel of peak emission.
;
; CATEGORY:
;       distance-omnibus
;
; CALLING SEQUENCE:
;       BGPS_POSTAGE [,BGPS_ML=bgps_ml][,RADIUS=radius][,/PLOT]
;                    [,START=start][,REAR=rear]
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       BGPS_ML -- IDL Save file with BGPS Map Locations [Default:
;                  local/bgps_v102_map_locations.sav]
;       RADIUS  -- Radius of the postage stamp image (in pixels).
;                  [Default: 35 -- final image is 8.52' x 8.52']
;
;       *Note: (CNUM_LIST) and (START, REAR) are mutually exclusive
;              options: if CNUM_LIST is supplied, then START & REAR
;              are ignored. 
;       CNUM_LIST -- List of BGPS catalog numbers for objects to
;                    validate (as opposed to running the entire
;                    catalog).  (Also, sets /PLOT keyword.)
;       START     -- First BGPS catalog number to process
;                    [Default: #1]
;       REAR      -- Last BGPS catalog number to process [Default:
;                    last entry in the catalog.
;
; KEYWORD PARAMETERS:
;       /PLOT   -- Display plots to screen rather than writing .FITS files.
;
; OUTPUTS:
;       Writes FITS file of each BGPS source to the directory
;       ./local/postage  If keyword /PLOT is set (i.e. plot to
;       screen), then FITS files will not be generated.
;
; OPTIONAL OUTPUTS:
;       NONE

; COMMON BLOCKS:
;       NONE
;
; EXAMPLE:
;       bgps_postage,radius=35,start=8000,/PLOT
;       bgps_postage
;
; MODIFICATION HISTORY:
;       Created:  08/27/10, TPEB -- Finally codified and documented
;                                   this procedure.
;       Modified: 04/29/11, TPEB -- Changed keyword NOPLOT -> PLOT, so
;                                   that the default action is to
;                                   create the .FITS files.  Also,
;                                   added HISTORY to FITS header
;                                   including routine name and date of
;                                   creation.  Plus, code cleanup.
;       Modified: 07/18/11, TPEB -- Added FITS keyword ORIGFILE and
;                                   FILENAME to postage-stamp files to
;                                   make finding the originating file
;                                   easier.
;       Modified: 07/19/11, TPEB -- Added check for correct PPBEAM
;                                   FITS keyword in the BGPS header.
;       Modified: 09/12/11, TPEB -- Updated START/REAR/CNUM_LIST to
;                                   conform in usage with other
;                                   routines in the distance-omnibus
;                                   project.
;       Modified: 12/09/11, TPEB -- Modified way postage stamp
;                                   subsection is defined (also
;                                   removes conflict with keyword
;                                   START).
;
;-

PRO BGPS_POSTAGE, BGPS_ML=bgps_ml, RADIUS=radius, PLOT=plot, $
                  START=start, REAR=rear, CNUM_LIST = cnum_list
  
  COMPILE_OPT OBSOLETE
  
  ;; Parse keyword parameters
  IF ~ KEYWORD_SET( radius ) THEN radius = 35
  do_plot = KEYWORD_SET( plot )
  IF ~ KEYWORD_SET( bgps_ml ) THEN $
     restore,'local/bgps_v102_map_locations.sav',/verbose ELSE $
        restore,bgps_ml,/ver
  
  ;; Set up the plot environment, if desired
  IF do_plot THEN BEGIN
     set_plot,'x'
     window,2
     cgLoadct,3
  ENDIF
  
  ;; Read in the BGPS Distance Database CSV file
  s = read_bgps_csv('bgps_distance_database.csv',csv,/VERBOSE)
  
  ;; If CNUM_LIST supplied, then only do the big loop for those objects
  IF n_elements( cnum_list ) NE 0 THEN BEGIN
     n_bgps = n_elements( cnum_list ) 
     ind = WHERE_ARRAY(cnum_list, s.cnum)
     s = s[ind]
     start = 0L
     rear = n_bgps-1
  ENDIF ELSE BEGIN
     ;; Else, look for START & REAR keywords and adjust accordingly
     n_bgps = csv.nrows
     IF n_elements(start) NE 0 THEN start = long(start) ELSE start = 0L
     IF n_elements(rear) NE 0 THEN rear = long(rear) ELSE rear = n_bgps-1
  ENDELSE  
  
  ;; Set up directory locations, and get BGPS Map and Label files
  LABEL_DIR  = 'local/bgps/label/'
  SMOOTH_DIR = 'local/bgps/other/'
  MAP_DIR    = 'local/bgps/maps/'
  POSTAGE_D  = 'local/postage/'
  bgps_list  = FILE_SEARCH(LABEL_DIR +'*.fits.gz', COUNT=n_list)
  bgps_maps  = FILE_SEARCH(MAP_DIR   +'*.fits', COUNT=n_map)
  bgps_smoo  = FILE_SEARCH(SMOOTH_DIR+'*_13pca_smoothmap50.fits', COUNT=n_smoo)
    
  q = '"'
  
  print,n_list,n_map,n_smoo
  IF (n_list NE n_map) OR (n_smoo NE n_map) THEN $
     message,'Mismatch between maps and labels!'
  
  ;; Read the map images & headers to memory (for faster processing)
  FOR j=0L, n_list-1 DO BEGIN
     
     command = string(j,q,bgps_maps[j],q,j,format=$
                      "('image',I0,' = readfits(',A0,A0,A0,',hdr',I0,')')")
     errcode = Execute(command)
     
     command = string(j,q,bgps_list[j],q,j,format=$
                      "('label',I0,' = readfits(',A0,A0,A0,',lhdr',I0,')')")
     errcode = Execute(command)
     
     command = string(j,q,bgps_smoo[j],q,j,format=$
                      "('smooth',I0,' = readfits(',A0,A0,A0,',shdr',I0,')')")
     errcode = Execute(command)
     
  ENDFOR
  
  plotj = -1
  
  ;;Loop through BGPS sources
  FOR i=start, rear DO BEGIN
     
     ;; Get image index # for extraction
     j = WHERE(bgps_maps EQ MAP_DIR+bgps[i].filename, nj)
     
     ;; Speed up processing by using settings from previous loop, if possible
     IF (j NE plotj) THEN BEGIN
        command = string(j,format="('image = image',I0)")
        errcode = Execute(command)
        
        command = string(j,format="('label = label',I0)")
        errcode = Execute(command)
        
        command = string(j,format="('smooth = smooth',I0)")
        errcode = Execute(command)
        
        command = string(j,format="('hdr = hdr',I0)")
        errcode = Execute(command)
     
        command = string(j,format="('lhdr = lhdr',I0)")
        errcode = Execute(command)

        command = string(j,format="('shdr = shdr',I0)")
        errcode = Execute(command)
        
        ;; Also, check the FITS header for the correct value of PPBEAM,
        ;;   and write the corrected value, if necessary.
        IF check_bgps_ppbeam( hdr ) THEN writefits,bgps_maps[j],image,hdr
        
        IF do_plot THEN BEGIN
           plotrange = set_plot_range(image)
           plotimage,image,range=plotrange,/preserve_aspect,$
                     title=bgps[i].filename
        ENDIF
        
        plotj = j
     ENDIF
     
     ;; Create Plot title, and subsection parameters for this BGPS source
     titlestr = string(s[i].cnum,s[i].glon_peak,s[i].glat_peak,format=$
                       "('BGPS #',I4,'  l = ',F5.1,'  b = ',F5.2)")
     
     ;; Determine box corners for postage stamp extraction
     imgsz = size(image,/DIM)
     box = long( [ (bgps[i].xpos-radius)>0, $
                   (bgps[i].xpos+radius)<(imgsz[0]-1), $
                   (bgps[i].ypos-radius)>0, $
                   (bgps[i].ypos+radius)<(imgsz[1]-1) ] )
     
     ;; Use HEXTRACT to extract the subsection, and update header WCS
     hextract,image,hdr,ps,pshdr,box[0],box[1],box[2],box[3],/silent
     hextract,smooth,shdr,pss,psshdr,box[0],box[1],box[2],box[3],/silent
     hextract,label,lhdr,pl,plhdr,box[0],box[1],box[2],box[3],/silent
     
     ;; Make the pl image just the source mask (1 - 0)
     pl = pl EQ bgps[i].cnum
     
     ;; Write the postage stamp (image and label) to FITS files  (but
     ;; only if not plotting to screen)
     IF ~ do_plot THEN BEGIN
        
        ;; Add BGPS CNUM & Object Name and HISTORY to the FITS headers
        ;; including date of modification and the routine name.
        
        hist = ['BGPS_POSTAGE: '+systime(0),$
                'Created BGPS postage-stamp image']
        GET_DATE,dte
        
        ;;==============================================================
        ;; DO 3 times, on for each postage-stamp image created...
        psfn  = 'bgps'+string(s[i].cnum,format="(I04)")+'.fits'
        pssfn = 'smooth'+string(s[i].cnum,format="(I04)")+'.fits'
        plfn  = 'label'+string(s[i].cnum,format="(I04)")+'.fits'
        
        
        ;; Smoothed image
        sxaddpar, psshdr, 'DATE_PS',dte, $
                  ' Creation UTC date of postage-stamp FITS file', AFTER='DATE'
        sxaddpar, psshdr, 'OBJNAME' ,s[i].name, $
                  ' Object Name in Galactic Coordinates', AFTER='DATE_PS'
        sxaddpar, psshdr, 'CNUM', s[i].cnum, $
                  ' Bolocat (v1.0.1) catalog number', AFTER='OBJNAME'
        fn = strmid(bgps_smoo[j],strpos(bgps_smoo[j],'/',/REVERSE_SEARCH)+1)
        sxaddpar, psshdr, 'ORIGFILE', fn[0], $
                  ' Original BGPS file', AFTER='CNUM'
        sxaddpar, psshdr, 'FILENAME', pssfn, AFTER='ORIGFILE'
        sxaddhist, hist, psshdr
        
        ;; BGPS science image
        sxaddpar, pshdr, 'DATE_PS',dte, $
                  ' Creation UTC date of postage-stamp FITS file', AFTER='DATE'
        sxaddpar, pshdr, 'OBJNAME' ,s[i].name, $
                  ' BGPS Object Name in Galactic Coordinates', AFTER='DATE_PS'
        sxaddpar, pshdr, 'CNUM', s[i].cnum, $
                  ' Bolocat (v1.0.1) catalog number', AFTER='OBJNAME'
        fn = strmid(bgps_maps[j],strpos(bgps_maps[j],'/',/REVERSE_SEARCH)+1)
        sxaddpar, pshdr, 'ORIGFILE', fn[0], $
                  ' Original BGPS file', AFTER='CNUM'
        sxaddpar, pshdr, 'FILENAME', psfn, AFTER='ORIGFILE'
        sxaddhist, hist, pshdr
        
        ;; BGPS label map
        sxaddpar, plhdr, 'DATE_PS',dte, $
                  ' Creation UTC date of postage-stamp FITS file', AFTER='DATE'
        sxaddpar, plhdr, 'OBJNAME' ,s[i].name, $
                  ' Object Name in Galactic Coordinates', AFTER='DATE_PS'
        sxaddpar, plhdr, 'CNUM', s[i].cnum, $
                  ' Bolocat (v1.0.1) catalog number', AFTER='OBJNAME'
        fn = strmid(bgps_list[j],strpos(bgps_list[j],'/',/REVERSE_SEARCH)+1)
        sxaddpar, plhdr, 'ORIGFILE', fn[0], $
                  ' Original BGPS file', AFTER='CNUM'
        sxaddpar, plhdr, 'FILENAME', plfn, AFTER='ORIGFILE'
        sxaddhist, hist, plhdr
        
        ;;==================================
        ;; Write the FITS files
        mwrfits,pss,POSTAGE_D+pssfn,psshdr,/create
        mwrfits,ps,POSTAGE_D+psfn,pshdr,/create
        mwrfits,pl,POSTAGE_D+plfn,plhdr,/create ;,/silent
        
     ENDIF
     
     ;; Plot to screen, if desired.  Includes wait time for visual inspection
     IF do_plot THEN BEGIN
        plotimage,ps,range=set_plot_range(ps*(pl*0.65 + 0.35)),$
                  /preserve_aspect,/interp,$
                  title=titlestr,xstyle=4,ystyle=4
        
        ;; Plot contours
        xc = findgen(subsize[0]) + 0.5
        yc = findgen(subsize[1]) + 0.5
        contour,pl,xc,yc,levels=[0.5],thick=3,color='ffff00'x,/overplot
        plots,radius+0.5,radius+0.5,psym=2,thick=3,symsize=2,color='22DD22'x
        
        ;; Plot coordinate axes
        plot_wcs_axes, pshdr, 1
        
        wait,0.5
     ENDIF
     
  ENDFOR
  
END
